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[1] The paper is devoted to the development of the asymptotic algorithm for the cloud top 
height A and the geometrical thickness / determination using measurements of the cloud 
reflection function. It is based on the asymptotic theory of the radiative transfer in the 
oxygen absorption bands and simple parameterization of the radiative transport in the 
atmosphere above and under a cloud. In particular, we have studied the influence of 
the error of the developed approximate theory on the accuracy of the retrieval of the pair 
(h, L). It was assumed that there is only a single cloud layer having the same value of the 
liquid water content in all points inside the cloud. The values (h, /) have been found, 
solving the inverse problem having as input the synthetic spectra of backscattered light. 
The synthetic spectra were found using the exact solution of the forward problem for given 
values of (A, L). The retrieval technique was based on the asymptotic theory. We have 
found that the error of the cloud top height determination is smaller than 20 m, and the 
error of the cloud geometrical thickness determination is smaller than 500 m for solar 
angles 20—70 degrees, values of A in the range 1—12 km, values of geometrical thickness / 
in the range 0.5—2 km, and values of the cloud optical thickness changing in the range 
10-50. The surface albedo has been assumed to be equal zero. We have also studied the 
influence of the cloud liquid water profile on the results of the retrieval of the pair (A, L). It 
was found that the error of the cloud top height determination increases up to 600 m, 

if the assumed cloud has a changing with height liquid water content, and retrievals are 
made applying the inversion with assumed constant liquid water content profile. The error 
in the geometrical thickness increases up to | km in this case. Errors of the retrieval 
increase even further if the retrieval, on the basis of the homogeneous layer theory, is 
applied to the two-layered cloud system (e.g., the upper cloud consists of ice crystals and 
the lower cloud is in a liquid phase). This signifies the importance of the cloud vertical 
inhomogeneity on the retrieval results. | INDEX TERMS: 0305 Atmospheric Composition and 
Structure: Aerosols and particles (0345, 4801); 0320 Atmospheric Composition and Structure: Cloud physics 
and chemistry; 0343 Atmospheric Composition and Structure: Planetary atmospheres (5405, 5407, 5409, 
5704, 5705, 5707); 0345 Atmospheric Composition and Structure: Pollution—urban and regional (0305); 
KEYWORDS: radiative transfer, clouds, remote sensing 
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1. Introduction 


[2] Cloud altitude and geometrical thickness are impor- 
tant parameters for a number of meteorological and clima- 
tological applications. For instance, the very position of a 
cloud top height may be an indication of an inversion layer 
in the atmosphere. Cloud altitudes and types indicate the 
thermodynamic and hydrodynamic structure of the atmo- 
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sphere. These parameters affect energy budgets and radia- 
tive heating. In particular, higher altitudes of clouds lead to 
higher surface temperatures. 

[3] The distribution of heating in the cloud, which is 
influenced by the cloud top height and the geometrical 
thickness, is of importance for the cloud dynamics and the 
evolution of cloud microstructure. Also, cloud top and 
cloud base heights influence the photon horizontal transport 
in clouds [Titov, 1998]. This has an importance to a number 
of issues and in particular to a so-called cloud absorption 
anomaly problem [Stephens and Tsay, 1990]. 

[4] This emphasizes the importance of information on the 
cloud top height A, the cloud geometrical thickness /, and 
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the cloud base height b = h — / derived on a global scale. 
Values h and / can be, in principle, retrieved using satellite 
measurements. A great number of methods has been pro- 
posed for the determination of the cloud top height A. They 
range from active measurements [see, e.g., Winker et al., 
1999; Poole et al., 2003], using lidars and radars on satellite 
platforms, to passive techniques, on the basis of the pro- 
cessing of information, contained in the thermal infrared 
radiances [King et al., 1992], ring effect [Joiner and 
Bhartia, 1995; de Beek et al., 2001], reflected light polar- 
ization [Knibbe et al., 2000], and stereo photogrammetric 
technique [Moroney et al., 2002], to name a few. 

[s] In particular, measured thermal infrared radiances are 
converted to cloud top pressures using assumed (or known 
from other sources) an atmospheric temperature profile 
[King et al., 1992; Platnick et al., 2003]. A possible error 
in this profile influences the retrieval results considerably 
[Naud et al., 2002]. 

[6] Most accurate estimations (the accuracy in A is better 
than 20 m) can be obtained from satellite laser systems. 
However, they are rather expensive and have a poor spatial 
coverage. 

[7] A great number of passive techniques developed also 
indicate that none of them can provide an accuracy compa- 
rable with laser active systems. In particular, Naud et al. 
[2002] found that Moderate resolution Imaging Spectroradi- 
ometer (MODIS), developed by NASA [King et al., 1992], 
gives the values of A, which differ from those obtained from 
the ground-based microwave measurements up to 1—2 km. 
MODIS data for h can differ from correspondent data 
obtained from the Multi-Angle Imaging Spectroradiometer 
(MISR) by 2 km. Note that MODIS low cloud top height 
product is based on the brightness temperature measured at 
11 um. For high clouds, MODIS algorithm is relied upon the 
CO>-slicing method [King et al., 1992]. Then the CO2 15 um 
absorption band is used for the retrieval procedure. The 
MISR cloud altitude retrievals are derived from multiangle 
red channel radiances and a stereo photogrammetric tech- 
nique [Moroney et al., 2002]. 

[s] Generally, measurements in visible and near-infrared 
cannot be used to retrieve the parameters A and /. This is due 
to an extremely weak sensitivity of the top-of-atmosphere 
radiance / to these parameters in this spectral region. Indeed, 
the radiative transfer equation, which is usually applied to the 
interpretation of satellite measurements, has as a parameter 
the value of the optical thickness T, which is a product of / and 
a priori unknown value of the cloud extinction coefficient ext 
(for homogeneous clouds). The value of t can be easily 
obtained from measurements in visible [Nakajima and King, 
1990] or near-infrared [Platnick et al., 2003]. Measurements 
in near-infrared can be also used to find the effective radius of 
particles aef [Nakajima and King, 1990; Nakajima et al., 
1991; Kokhanovsky et al., 2003] and cloud thermodynamic 
state [Knap et al., 2002]. However, the information on the 
parameters / and / is virtually absent. This is actually of a 
benefit for the retrievals of dep , T, and liquid water path W, 
which is proportional to the product of agp and T. Note that it 
follows approximately: 


W= AdefT, 


where A = 2p/3 and p = 1 gem ° is the density of water. 
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[9] The situation is radically changed if we consider the 
radiative transfer in the molecular absorption line [ Yamomoto 
and Wark, 1961; Heidinger, 1998; Stephens and Heidinger, 
2000; Heidinger and Stephens, 2000, 2002]. Indeed, let us 
assume that we have a gas in a planetary atmosphere, which 
absorbs almost all incident radiation in a narrow band. Then 
the depth of this band, measured by a receiver on a satellite, 
will depend on the cloud altitude. Gas concentrations gener- 
ally decrease with the distance from the ground. Therefore 
clouds at a high altitude do not allow photons to penetrate to 
low atmospheric layers and be absorbed there. So the depth of 
a molecular line will decrease if high clouds are present in the 
field of view of a sensor. 

[10] This idea for the cloud top height retrieval was 
proposed by Hanel [1961] and has already been applied 
for cloud top altitude measurements. In particular, Saiedy et 
al. [1967] analyzed measurements obtained by the Gemini5 
astronauts in the region of the oxygen A band centered at 
760 nm. Measurements were performed using a small hand- 
held spectrograph with the spectral resolution 5 A or 10 A, 
depending on settings. They were able to determine altitude 
of several cloud systems (e.g., Sc cloud over the east Pacific, 
a cloud in the Intertropical Convergence Zone, the hurricane 
Doreen, etc.) [see Saiedy et al., 1967, Table 1]. In particular, 
the absorption line for the area over a hurricane was not so 
pronounced as for a case of a low Sc cloud. They found that 
the hurricane top pressure p was 323 mb (A ~ 7 km). The 
value of p for a Sc cloud was 835 mb (h ~ 2 km). 

[11] An important finding of their work was the fact that 
the cloud geometrical thickness / cannot be ignored in the 
retrieval procedure. Saiedy et al. [1967] emphasized that the 
measured cloud top altitudes, on the basis of the assumption 
that the cloud is substituted by the Lambertian diffuser with 
zero transmittance, give values that are always below actual 
cloud top heights. 

[12] This can be easily understood on physical grounds as 
well. Indeed, in reality photons penetrate through a cloud 
top in inner cloud areas (and also below a cloud) and are 
absorbed by oxygen there. This leads to the increase of the 
oxygen absorption depth even for high clouds. If the 
process of the photon penetration is neglected (a Lambertian 
diffuser assumption), then this increased depth of the 
absorption line is interpreted as an existence of a cloud on 
a level that is lower that the actual altitude of a cloud. This 
actually was observed by Saiedy et al. [1965] while com- 
paring retrieved cloud top heights (from airborne reflectance 
spectra) with those obtained using airborne in situ measure- 
ments of h. They also introduced a correction factor 
(typically 1.06—1.38) to account for this effect. Their 
findings lead to several conclusions. 

[13] First of all, for a correct cloud top height retrieval, 
one should fully account for a photon transport and multiple 
scattering inside a cloud. This is confirmed also by a recent 
study of Vanbauce et al. [2003]. Therefore the substitution 
of a real cloud by a Lambertian reflector model leads to 
potentially large retrieval errors, depending on an actual 
cloud geometrical thickness. 

[14] Second, the cloud geometrical thickness / can be 
estimated from independent satellite measurements [Asano 
et al., 1995; Kuji and Nakajima, 2002]. Then retrieved 
cloud geometrical thickness might be used to obtain a 
realistic estimation of the light absorption inside a cloud 
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and therefore to reduce uncertainties in the cloud top 
altitude measurements. Such a suggestion was made for 
the first time by Saiedy et al. [1965]. Following this path, 
Asano et al. [1995] and Kuji and Nakajima [2002] used the 
oxygen A band for the cloud geometrical thickness deter- 
mination for several case studies. Note that Hayasaka et al. 
[1995] used the solar reflection at 0.94 um water vapor 
absorption band to estimate /. 

[15] The development of the radiative transfer theory, an 
accurate oxygen absorption vertical profiling, and rapid 
progress in computer technology allowed to solve the 
problem of cloud top height determination using O2 A band 
for a known value of the cloud geometrical thickness (at 
least for homogeneous cloud layers) [see, e.g., Fischer and 
Grassl, 1991]. For this, one should run the radiative transfer 
code for underlying surface-cloudy atmosphere system for 
lines inside the oxygen A band and minimize the difference 
between measured and simulated oxygen absorption band 
spectra. The cloud optical thickness, cloud thermodynamic 
state, and effective radius of particles, which are also needed 
in fitting procedure, can be obtained from measurements 
outside gaseous absorption bands as specified, e.g., by 
Nakajima and King [1990] and Knap et al. [2002]. In a 
similar way the cloud geometrical thickness can be estimated 
(for a known cloud top height, e.g., from space lidar 
measurements). 

[16] The inverse problem can be also posed for a simul- 
taneous determination of cloud top height and cloud geo- 
metrical thickness. For this, however, a high spectral 
resolution and multispectral measurements in the oxygen 
A band are needed. The alternative is the use of the oxygen B 
band centered at 687 nm, other molecular absorption bands, 
or their combinations. 

[17] It should be stressed that exact radiative transfer 
calculations cannot be applied for satellite operational 
algorithms of the cloud top height determination owing to 
high requirements in respect to a computational speed. This 
prompts for a use of a look-up table approach. However, the 
estimation of the size of such a look-up table gives at least 
10° cases. The size of the database can be in principle 
reduced using the polynomial approach [Fischer and 
Grassl, 1991] or neural networks techniques [Fischer et 
al., 2000]. 

[1s] The main shortcoming of these methods is, however, 
inflexibility. For instance, the change of the position, width, 
or the number of spectral channels forces to construct a new 
database or neural network training. This does not allow to 
use algorithms developed for a given spectrometer or 
radiometer to interpret data obtained from different instru- 
ments on board multiple satellite platforms using the same 
database constructed. Similar problems arise also if one 
needs to feed the algorithm with updated spectroscopic 
information. 

[19] A different approach to reducing the size of the 
database is based on the use of the asymptotic radiative 
transfer theory [van de Hulst, 1980]. This allows at the same 
time to preserve a wide range of the model applicability in 
respect to cloud parameters and illumination/observation 
geometry. This is similar to numerical methods of the 
radiative transfer equation solution in this respect. Such a 
technique is used, e.g., in the MODIS cloud retrievals for 
the estimation of cloud optical thickness and the droplet/ 
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crystal size determination [King et al., 1997] if the cloud 
optical thickness is larger than 10. Then the accuracy of the 
asymptotic theory is better than 1%, which is smaller then 
possible error in a forward model. 

[20] We propose to extend this technique of the database 
reduction to the molecular absorption bands. To this end, 
modified asymptotic equations [Kokhanovsky and Rozanov, 
2003, 2004], which have an accuracy typically better than 
5% for a cloud optical thickness larger than 5, can be used. 
The exploration of such a possibility for the problem of the 
retrieval of cloud geometrical parameters (A, /) is a main 
task of this paper. 

[21] Thinner clouds can be in principle handled by a look- 
up table approach. However, such cloud systems are highly 
inhomogeneous in a horizontal direction, making the appli- 
cation of the standard one-dimensional radiative transfer 
equation not justified. Therefore we consider only the case 
of optically thick clouds here. 


2. Forward Model 


2.1. Radiative Transfer Model: 
The Monochromatic Case 


[22] The analytical radiative transfer forward model, 
which is used here to develop the cloud top height and 
geometrical thickness retrieval algorithm, is described by 
Kokhanovsky and Rozanov [2003]. They also studied the 
accuracy of the model, which is typically better than 5% for 
the calculation of the reflection function at the top of 
atmosphere (at the cloud optical thickness larger than 5). 

[23] For the sake of completeness we formulate the model 
here as well without any derivations, which can be found 
elsewhere [Kokhanovsky and Rozanov, 2003]. The geome- 
try of the problem is given in Figure 1. Unpolarized solar 
light illuminates the aerosol-gaseous cloudy atmosphere in 
the direction Qo (Uo, Yo), and a receiver detects the reflected 
light intensity / in the direction, specified by Q (9, p). Here 
Ùo = arccos € and 0 = arccos |n| are zenith incident and 
observation angles, respectively; Yo and » (not shown in 
Figure 1) are correspondent azimuth angles. We will assume 
that po = 0. Then the azimuth difference between incident 
and reflected beams directions is given by the angle y. 

[24] A single horizontally homogeneous plane-parallel 
cloud layer is characterized by a top height A and the 
geometrical thickness / (see Figure 1). Both 4 and / can 
be in principle arbitrary. However, we assume that / > 52 (or 
T > 5), where £ = os. is the photon free path length in a 
cloud. For a typical value of £ = 20 m, it follows: Z > 100 m. 
This gives an approximate lower value of the cloud geomet- 
rical thickness, which can be handled in the framework of 
our model. Cloud top height can vary from approximately 
0.1 km till the top of atmosphere (TOA) in the model. 
However, note that terrestrial water and ice clouds, having 
large optical thickness, do not penetrate to heights larger 
than approximately 10—15 km, depending on the latitude. 

[25] Vertical profiles of the gaseous and aerosol absorp- 
tion and scattering are fully accounted for in the model as 
described by Kokhanovsky and Rozanov [2003]. We also 
account for a vertical variation of a cloud single scattering 
albedo, taking into consideration the increase of gaseous 
absorption toward the cloud base. The vertical variation of 
the cloud liquid content is fully accounted for as well. 
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Figure 1. 


However, the phase function is assumed to be the same in 
all parts of a cloud. It is found using the Mie theory 
[Kokhanovsky, 2001] for the effective radius (the ratio of 
the third to the second moment of the droplet size distribu- 
tion) equal to 6 um as in Deirmendjian’s [1969] cloud C1 
model. 

[26] The TOA monochromatic reflection function R(\, 
Ùo, Ù, p) is related to the reflected monochromatic light 
intensity J(\, Ùo, Ù, p) by the following relationship 
[Thomas and Stamnes, 1999]: 


os I(, 0, Ùo, 9) 


RON, 0, 00:9) = F(X) cos 9 ’ (o) 


where nF (A) cosùọ is the incident spectral solar flux 
density on the TOA. The value (1) can be easily obtained 
from measurements performed by various radiometers and 
spectrometers placed on satellite platforms [Liou, 1992]. 

[27] Following the results given by Kokhanovsky and 
Rozanov [2003], the TOA reflection function R is presented 
as a sum of two terms (see Figure 1): 


R=R,+T RT, (2) 


where R, describes light scattering and radiative transfer in 
the atmosphere above a cloud (with account for both 
gaseous and particle scattering and absorption). The value 
of R, is due to cloud-underlying atmosphere and surface 
contribution (see Figure 1). The multiplier 7; accounts for 
the extinction of the direct solar light on the path from the 
Sun to the cloud top, and T, accounts for the same effect 
but on the way from a cloud to a satellite receiver (see 
Figure 1). We omit arguments in functions in equation (2) 
for the sake of simplicity. 

[28] The scattering of light above clouds is rather weak. 
Thus the value of R, is calculated in the single scattering 
approximation (see Appendix A). 


The geometry of the problem. 


[29] It should be pointed out that coupling between 
atmospheric layers above and below cloud top is neglected 
in equation (2). To account partially for this coupling and 
multiple light scattering above the cloud top, we assume for 
a product T = 7,;7,: 


T= exp[—1/(€7! + a) (3) 


where 
N H 
T= > J Sre (4) 
h 


Here H = 60 km is the assumed TOA height, Chri) is the 
ith gas absorption cross section, N is the total number of 
gases present, c; is the ith gas concentration, and £ = cos Uo, 
n = |cos |. Therefore only extinction owing to gaseous 
absorption contributes in 7’ in the framework of our forward 
model. 

[30] Actually, the value of tr’ should include also extinc- 
tion by aerosol particles (the optical thickness T^) and 
molecules (the optical thickness 7®). We neglect them in 
equation (3), which results in increase of T (and 7;R,T,. in 
equation (2)). This partially compensates for multiple light 
scattering in the layer above cloud, which otherwise is 
completely neglected in this study. Note that this is a 
standard method for a partial account of multiple light 
scattering effects in the problems of the atmospheric cor- 
rection. However, the definition of 7’ in equation (3) differs 
depending on the problem and a spectral range studied [see, 
e.g., Gordon et al., 1983; Wrigley et al., 1992; Wang and 
King, 1997]. 

[31] Now we consider the below—cloud top contribution 
R,. The general expression for this term for the cloud- 
underlying Lambertian surface with the albedo A is well 
known [Liou, 1992]. It has the following form, if we neglect 
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the contribution of the direct light, which is a valid 
assumption for thick clouds considered here: 


Ate(9o)te(9) 


Ry (9, vo, p) _ RC, Vo, p) pa 
1 — Ar, 


(5) 


where Re (Vo, Ù, p) is the cloud reflection function at A = 0, 
te (0) is the diffused transmittance of a cloud layer, and re is 
the spherical albedo of a cloud. Convenient approximate 
equations for Re (Ùo, Ù, $), te (9), and r, are summarized in 
Appendix A. 

[32] Equations (2), (3), (5), and those given in Appendix A 
do not completely solve the problem we face. Indeed, there is 
aerosol-gaseous medium between underlying surface and a 
cloud base (see Figure 1). It will influence the relationship 
(5) obtained for the case of a transparent layer between a 
cloud and a ground surface. 

[33] To approximately account for this influence, we 
substitute A in equation (5) by the effective albedo 
[Kokhanovsky and Rozanov, 2003] 


At 


At = Na kj 
at Ar 


(6) 


where r, is the spherical albedo of the aerosol layer always 
existing between a cloud and an underlying surface (see 
Figure 1) and ¢, is the aerosol total transmittance. We will 
neglect the aerosol absorption at this point and assume that 
ta = 1 — ra. Thus for the calculation of A* we should know 
the value of r,. Accounting for the fact that the aerosol 
optical thickness is usually small (smaller than 0.1—0.3) and 
the aerosol contribution to the TOA reflectance from 
underneath of thick clouds is weak, we expect that rather 
coarse approximations for the calculation of r, can be used. 
In particular, we have calculated r, in the framework of the 
single scattering approximation, described by Wiscombe 
and Grams [1976] (see Appendix A). For thicker aerosol 
layers the approximation developed by Kokhanovsky and 
Mayer [2003] can be used. 

[34] To complete the model, we also account for the light 
absorption of diffused light fluxes by gases present in a layer 
between a cloud and an “effective” underlying surface. To 
do so, we multiply ¢ (o) and t- (Ù) in equation (5) by the 
flux transmittance [Feigelson, 1984; Li, 2000]: 


1 
y=2 as sop- (7) 
jasc 


where T* is the optical thickness due to absorbing gases in a 
layer between a cloud bottom and a ground surface. Note 
that it follows with a high accuracy [Kokhanovsky and 
Rozanov, 2003]: y © exp (—1.7t*) at T* < 0.4. 

[35] Summing up, we have the following final approxi- 
mate relationship for the TOA reflection function: 


A*y"t.(00)te(0) A 


RO, Do, p) = R (ð, Do, P) za R.(0, Bo, P) + 1 — A*r 


(8) 


where A* is given by equation (6) and approximate analytical 
formulae for Ra (Ù, Yo, P), Re (9, Yo, P), Fe te (Do) are 
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presented in Appendix A. In conclusion, we emphasize that 
our model accounts for (see Figure 1): (1) light extinction 
on paths from Sun to a cloud and from a cloud to a satellite; 
(2) multiple photon scattering inside a cloud; (3) single 
photon scattering above a cloud; (4) surface and aerosol light 
reflection and absorption beneath a cloud; and (5) vertical 
variation of scattering and absorption characteristics of gases 
and particulate matter both inside and outside ofa cloud. Note 
that some elements of the advanced forward model presented 
here can be used far beyond a rather narrow topic considered 
in this paper. 

2.2. Radiative Transfer in a Finite Spectral Interval: 
The Correlated k Distribution Scheme 


[36] The TOA intensity Za}, measured by a radiometer or 
spectrometer, is in fact an integral of the monochromatic 
intensity /(\, Ùo, 0, p) taken with respect to the wavelength. 
It can be presented by the following equation: 


M+L 


Im (ùo, 0,9) = J 
NL 


SONION, 90,9, 9)dN', (9) 


where A» is the spectral resolution of the instrument and 


fO, N) is the instrument spectral response function, 


normalized as 


+L 
J fON)dN =1. 


A-L 


(10) 


It is assumed that f(\, \’) differs from zero only in the 
interval X € [\ — L, \ + L], where X is the central 
wavelength of the instrument response function. A similar 
integration should be performed to find the value of TF’, 
cos Ùo (the incident solar flux density on the TOA in the 
spectral range [\ — L, \ + L]), which is needed to determine 
the polychromatic reflection function 


_ Ly (9, 90, p) 


Ra (9, 00, 9) = Fosh (11) 


[37] The shape and the half-width of the function fA, \’) 
is only of a minor importance for measurements outside of 
atmospheric gaseous absorption bands. This is not the case, 
of course, for absorption bands of atmospheric gases (e.g., 
Oz, H20, CO2, N20, CHy, etc.), which have an extremely 
complicated structure. The response function fA, X^) 
depends on the particular instrument. For the theoretical 
study, this function can be taken in the following form: 


1/AX, X € [A — A)/2, A + Ad/2] 
f(N,») = , (12) 


0, in other cases 


which is a simple box function. Then we have L = AW/2. 
Also, one can use the Gaussian function: 


SAN) = (13) 


1 
Jana exp ( 
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Figure 2. The top-of-atmosphere spectrum measured by SCIAMACHY (points) for a cloudy scene over 
the Alps on 8 August 2002. Band widths of MERIS and GLI (shaded strips) are also shown. The 
resolution of GOME (almost the same as in SCIAMACHY) is not shown. 


where it is assumed that L — oo. In practical calculations 
we will assume that L = 3A. 
[38] It follows from equation (9) that 


Is ( (90,9, ¢) = Yori fA, XM (Xi, 90,9, P), (14) 


where <; are the quadrature coefficients and » is the central 
wavelength. Owing to the oscillating character and the 
complex structure of the function / (X) in the oxygen band, 
the number M should be large (typically 1000 at L = 0.5 nm). 
This is impractical, however, because the calculation of / 
(As Ùo Ù, p) using integrodifferential radiative transfer 
equation is a complex task, which takes a considerable time. 

[39] A so-called correlated k distribution scheme allows to 
reduce considerably the number of terms in equation (14). 
This is achieved replacing equation (14) by the following 
approximate expression: 


In (9o, 9,9) =S (A, N) SQV, 80, 0,9), (15) 


where m < M and S (\,, Ùo, 9, 9) = S Ddal, 0, y). The 


theory behind such a replacement i$ called the correlated 
k distribution scheme and described elsewhere [see, e.g., 
Lacis and Oinas, 1991; Buchwitz, 2000; Buchwitz et al., 
2000]. We have used constants D, given by Buchwitz 
[2000]. The values of J, (o, Ù, p) are calculated using 
the monochromatic radiative transfer equation for five 
precalculated profiles of the molecular oxygen absorption 
cross section (or exponential sum fitting coefficients 
[Buchwitz et al., 2000]), depending on the temperature, 
pressure, and the actual value of ^. Values of the oxygen 
absorption cross section used are given by Buchwitz [2000]. 


The correspondent database has been integrated in the 
SCIATRAN radiative transfer package [Rozanov et al., 
2002] and has been used in this work. 

[40] Buchwitz et al. [2000] proposed to use values of d; 
distributed from 755 to 775 nm (see Figure 2) with the step 
6 = 0.05 nm. This allowed to have an error smaller than 2% 
[Buchwitz, 2000] as compared to line-by-line calculations. 
Clearly we have for the value of m:m = [A)/6] for a box 
function (12). We have approximately for a Gaussian func- 
tion (13): m = [6A/6]. Therefore m depends on the spectral 
resolution of the instrument. Characteristic examples of 
values A of selected satellite instruments are presented in 
Table 1. In particular, we give data for the Global Imager 
(GLI) [Nakajima et al., 1998], POLDER [Hagolle et al., 
1999], the Global Ozone Monitoring Experiment (GOME) 
imaging spectrometer [Burrows et al., 1999], the Scanning 
Imaging Absorption Spectrometer for Atmospheric Chartog- 
raphy (SCIAMACHY) [Bovensmann et al., 1999], and the 
Medium Resolution Imaging Spectrometer (MERIS) [Bezy et 
al., 2000]. Note that GOME, SCIAMACHY, and MERIS are 
operating in space at the moment. 

[41] Figure 2 shows the spectral resolution of selected 
spectrometers as related to the oxygen A bandwidth. We 
see that both GLI and MERIS have much larger bandwidths 
as compared to SCIAMACHY and GOME instruments. 
Clearly the actual bandwidth and number of spectral points, 
where measurements are performed, influences the accuracy 
of the cloud top height retrieval in a great extent. 

[42] Let us substitute equation (8) in equation (15). Then 
we have for the proposed quasi-analytical model of the 
radiative transfer in a finite spectral interval: 


Ray (90,9, 9) = (X, y) 90,89), (16) 
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Table 1. Characteristics of Selected Space Instruments, Related to Measurements of the Backscattered Light in the 


Oxygen A Band (7550-7750 A) 


Instrument Platform Year Spectral Interval/Wavelength, A Spectral Resolution, A Spatial Resolution, km? 
GOME ERS-2 1995 5760—7940 3.3 40 x 320 
SCIAMACHY ENVISAT 2002 6040-8050 4.8 30 x 60 
MERIS ENVISAT 2002 7600 25 0.3 x 0.3 or 1.1 x 1.1 
GLI ADEOS-II 2002 7630 80 1.0 x 1.0 
POLDER ADEOS-II 2002 7633 100 6.0 x 7.0 

7651 400 


where R (Aj, Vo, 0, Y) = S DaRa(Vo, 9, p). The values of 
Ra(Ùo, Ù, p) are calculated at given effective absorption 
cross sections as discussed by Buchwitz et al. [2000]. 

[43] The crucial point in calculations according to equa- 
tion (16) is an account for the vertical distribution of the 
oxygen concentration (see, e.g., equation (4)). Also, one 
should account for the variation of the oxygen effective 
absorption cross section Cabs,a; With temperature and pres- 
sure. This is done using the database of Cabs,oj, calculated by 
M. Buchwitz (personal communication, 2003) for the values 
of à; in equation (15) using the HITRAN_ 2000 database 
[Rothman et al., 2003]. 

[44] The accuracy of equation (16) was studied by 
Kokhanovsky and Rozanov [2003] for a realistic vertically 
inhomogeneous atmosphere. One example of the perfor- 
mance of this equation is shown in Figure 3 for the GOME 
resolution with the instrument response function according 
to equation (13). We see that errors for the case considered 
are very low even inside the minimum of the absorption 
band, where the main assumption behind the derivation of 
the result for Re in equation (8) (namely weak light 
absorption) [see, e.g., Kokhanovsky et al., 2003] is violated. 


This is due to a comparatively large contribution of the term 
R, (small values of 7) in this case. 


3. Inverse Problem 
3.1. One-Parameter Retrieval Algorithm 

[45] Let us consider the inverse problem using equations 
introduced above. We assume at first that all parameters of 
our problem except the cloud top height A are known. Then 
the value of A can be found from spectral measurements in 
the oxygen A band using the following semianalytical 
technique. 


[46] Fist of all, we present equation (16) in the form of a 
Taylor expansion: 


R = R(ho) + Valh — ho)’, 
i=l 


(17) 


where a; = RO (ho)/l!, and we have omitted for the sake 
of simplicity the index A and angular variables (see 
equation (16)). Here R” (ho) is the / derivative of R at the 
point ho. The next step is the linearization, which is a 


0.6 E7 r — —J r r =i 
Eo SSS x h = 1.5 km 5 sy% J 
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Figure 3. The spectral reflection function of a cloudy medium in the oxygen A band calculated using 
equation (16) and the exact radiative transfer code at the cloud optical thickness equal to 20 and solar 
angle equal to 60° and nadir observation for the cloud C1 model [Deirmendjian, 1969] and the same 
aerosol and gaseous profiles as those specified by Kokhanovsky and Rozanov [2003]. Calculations 
correspond to the GOME resolution. The black underlying surface is assumed. 


7 of 21 


D05202 ROZANOV AND KOKHANOVSKY: SEMIANALYTICAL CLOUD RETRIEVAL D05202 
0.08 rrr prr rrr 0.08 rrr frrr rrr r 0.08 rrr prr yrr r 
= 0.5km = 1km | = 2km 
0.06 | J 0.06 F K J 0.06 : J 
0.04 0.04 0.04 F 
2 
S 0.02 0.02 0.02} 
E 
0.00 H --- 0.00 j 0.00 F- 
-0.02 -0.02 -0.02 } 
=0;04 Liai he eei S004 hiai hiii h =0 04 liate leje ii 
755 760 765 770 775 755 760 765 770 7755 755 760 765 770 7755 
Wavelength, nm Wavelength, nm Wavelength, nm 
Figure 4. The spectral dependence of derivatives “ and aR calculated using exact radiative transfer code 
(crosses) and our approximation (lines) for the same conditions as in Figure 3 except the solar angle equal 
to 53 degrees and cloud top height equal to 3 km for / = 0.5, 1, and 2 km. 
standard technique in the inversion procedures [Rozanov et where || || means the norm in the Euclid space of the 


al., 1998; Rodgers, 2000]. We found that the function R (A) 
is close to a linear one in a broad interval of the argument 
change [Kokhanovsky and Rozanov, 2003]. Therefore we 
neglect nonlinear terms in equation (17). Then it follows: 


R = R(ho) + R'(ho)(h — ho), (18) 


where R’ = R, We assume that R is measured at several 
wavelengths in the oxygen A band. Then instead of the 
scalar quantity R, we can introduce the vector Rmes with 
components (R (X1), R A2), ... R (A,)). The same applies to 
other scalars in equation (17). 

[47] Therefore equation (18) can be written in the follow- 
ing vector form: 

y= üi (19) 

where y = RŘmes— R (ho), @ = R' (ho) and x = h — ho. Note 
that both measurement and model errors are contained in 
equation (19). The solution x of the inverse problem is 
obtained by the minimizing the following cost function 
[Rodgers, 2000]: 


© = |7 — all’, (20) 


correspondent dimension. It follows from equation (20) 
[Rodgers, 2000] that 


(21) 


where (4, y) means the scalar product and £ gives the value 
of x, where the form (20) takes a minimum. Therefore, 
knowing values of the measured spectral reflection function 
Rmes and also values of the calculated reflection function 
R and its derivative R’ at h = ho and several wavelengths, 
the value of the cloud top height can be found from 
equation (21) and equality: h = x + ho. The value of ho can 
be taken equal to 1.0 km, which is a typical value for low 
level clouds. The main assumption in our derivation that the 
dependence of R on A can be presented by a linear function 
on the interval x [Kokhanovsky and Rozanov, 2003]. Note 
that the same equations can be used for the cloud 
geometrical thickness determination (with the substitution 
of h by /, see equation above). 

[48] Above, either value of h (or /) is assumed to be 
unknown. Our algorithm, however, has also a capability to 
find these two parameters (and also some other auxiliary 
information) simultaneously. This requires the minimization 
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aR and Æ & calculated using exact radiative transfer code 


(crosses) and our approximation (lines) for the same conditions as in Figure 4 except the cloud 


geometrical thickness is equal to 0.5 km and h = 


of the modified cost function (see equation (20)): & = 
|I} — AŽ ||”. The elements of the matrix A are correspondent 
weighting functions [Rodgers, 2000]. The solution of the 
inverse problem is given by the vector parameter X. The 
first two components of this vector give corrections to 
the initially assumed cloud top height and cloud top 
geometrical thickness. Others include information on aux- 
iliary parameters like the correction due to the assumed half- 
width of the spectrometer spectral response function. 

[49] Values of R and R’ are usually found from the 
numerical solution of the integrodifferential radiative trans- 
fer equation for a model atmosphere. They also can be 
calculated using equation (16). The last possibility allows to 
speed up the numerical solution of the inversion problem 
and will be used throughout the paper. 

[50] Note that it is possible to calculate derivatives R’ 
analytically. Correspondent expressions for R’ are cumber- 
some. They are not shown here. However, they were 
derived and used in the retrieval procedure. 

[sı] We compare numerically obtained from the exact 
radiative transfer code SCIATRAN [Rozanov et al., 2002] 
and approximate derivatives in Figures 4 and 5 It follows 
that derivatives (both with respect to h and /) differ from 
zero in the oxygen band. This actually allows for the values 
of h and / determination from the TOA reflectance measure- 


1, 3, 10 km. 


ments. It follows that analytical derivatives have a high 
accuracy. Their calculation is much faster as compared with 
the numerical perturbation technique, so we prefer their use 
in the retrieval procedure. 

[52] Note that Ry = = Ris positive and Rj = Ë is negative. It 
confirms that the 1 increase in the cloud top height leads to 
the increase of the reflectance at the fixed cloud geometrical 
thickness. The opposite is true for the increase of the cloud 
geometrical thickness at the fixed cloud top height. 

[53] An important step in every inverse problem is to 
understand the sensitivity of the model to the various 
possible errors [Rodgers, 2000]. They were studied in great 
detail by several authors using the exact radiative transfer 
theory [see, e.g., Wu, 1985; Fischer and Grassl, 1991; 
Fischer et al., 1991; O’Brien and Mitchell, 1992], so we 
do not repeat these studies here. In particular, Fischer et al. 
[1991] showed that all errors totaled up in the possible error 
in the cloud top height determination equal to 150 m as 
compared to in situ cloud top height measurements. The 
averaging procedure allowed then to reduce the error over 
the flight distance 20 km down to 40 m [Fisher et al., 1991]. 

[54] Clearly, if we use synthetic data obtained from the 
exact radiative transfer theory and also exact calculations 
for finding R and R’, then the error of finding h from 
equation (21) is close to zero. It is interesting to study how 
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Figure 6. The cloud top height error as the function of the actual cloud top height, obtained from the 
use of the analytical model presented here for known values of the cloud optical and geometrical 
thickness. Results for resolutions correspondent to MERIS and GLI instruments are shown (see text). The 
measured reflection function was simulated using the radiative transfer code SCIATRAN [Rozanov et al., 
2002] at the solar zenith angle equal to 60 degrees and nadir observations. 


large the error of the cloud top height determination is if we 
model measured data by the exact theory but retrieve the 
value of h, using approximate solutions for the reflection 
function and its derivatives introduced above. This shows us 
how errors of our forward model transform, e.g., in the error 
of the cloud top height determination A = h — h,, where h, is 
the retrieved value of the cloud top height. 

[s5] The dependence of A on the cloud top height for 
different values of /, T is given in Figure 6. The spectral 
resolutions chosen are correspondent to MERIS and GLI 
instruments (see Table 1). Note that for these instruments 
only a single-wavelength measurement in the oxygen A 
band is performed (see Figure 2). Then we have instead of 
equation (21): 


(22) 


where A, gives the retrieved value of the cloud top height. 
We use instead of measured value of the reflection function 
Rnes in equation (22), the correspondent function, which is 
obtained from the exact SCIATRAN [Rozanov et al., 2002] 
radiative transfer calculations. The measurement errors are 
neglected in this case. 

[s6] Note that for all retrievals performed in this paper, 
we have actually used values of Rmes, R, and R’ normalized 
to the average value of these functions outside the band (in 
the interval 755—757 nm (see Figure 3)). This enhances the 
accuracy of the retrieval. 

[57] It follows from Figure 6 that errors A of our 
technique for both instruments are smaller than 200 m, 
which is much better than the accuracy and which can be 
obtained from infrared passive techniques. Results for 
MERIS are slightly better because our approximation works 
better for smaller absorption [Kokhanovsky et al., 2003], 
which is the case for the MERIS measurements (see 
Figure 2). Note that for most of cases the value of A is 
negative, which means that the value of / is overestimated 


by the our approximate radiative transfer method. Therefore 
too high values of A are retrieved. 

[ss] Figure 7 gives the retrieved value of the cloud top 
height error A for different solar zenith angles. We see that 
the error is not influenced by the illumination condition 
dramatically. However, generally the error increases with 
the solar zenith angle. The influence of the solar zenith 
angle is more important for low clouds, where our technique 
gives larger retrieval errors. 

[s9] Although the error of the asymptotic algorithm (up to 
200 m) may appear unacceptably large (especially if other 
possible error sources are accounted for), it is small in 
comparison with common errors of passive methods used 
today for the cloud top height determination [Naud et al., 
2002]. They typically give the accuracy not better than 1km 
in the value of h [Naud et al., 2002]. 

[co] The question arises if the error of the retrieval can be 
reduced if the exact theory is applied to finding the cloud 
altitude from measurements performed in terrestrial atmo- 
sphere. As we stated above, the error A for synthetic data 
will be close to zero in this case. However, we believe that 
our analytical asymptotic cloud retrieval algorithm should 
give errors of A close to that obtained if the exact theory is 
applied to the interpretation of the experimental data 
obtained from measurements in the terrestrial atmosphere. 
This is mostly due to the fact that here plays a major role not 
the difference (approximately 5% or even less; see Figure 3) 
between the exact and asymptotic theory, but the uncertainty 
related to the forward model, which is basically the same 
both for the exact theory and approximation and does not 
necessarily correspond to the case presented during mea- 
surements (e.g., cloud horizontal inhomogeneity and the 
existence of several cloud layers). 

[61] However, this question can be cleared out only after 
such experiments are performed in situ, which is a complex 
task. We hope to make such a study in future. For this, 
collocated airborne and satellite measurements over extended 
cloud fields are clearly needed. As it was stated above, the 
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Figure 7. The same as in Figure 6 but for various solar zenith angles (SZA) at / = 0.5 km and T = 30. 


exact radiative transfer theory cannot be applied to opera- 
tional retrievals owing to speed requirements. 

[62] Some estimates of errors in the parameter h determi- 
nation for vertically inhomogeneous and realistic two-lay- 
ered cloud systems are given in section 4. To make them, 
we have performed a number of numerical experiments. We 
emphasize once more that the same algorithm as described 
above (after a slight modification (change A to /)) can be 
also used for the cloud geometrical thickness determination 
(e.g., if h is known from space lidar measurements). 


3.2. Two-Parameter Retrieval Algorithm 

[63] Generally, both MERIS and GLI (see Table 1) do 
not allow to find both parameters (h and /) from measure- 
ments in the oxygen A band only (one measurement). 
However, it is possible to obtain the cloud top height 
from the GLI infrared measurements [Nakajima et al., 
1998; Nakajima, 2001]. Then this information can be 
applied (at least in principle) to the determination of the 
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cloud geometrical thickness using the measurement in 
oxygen A band [Kuji and Nakajima, 2002]. However, a 
possible large error [Naud et al., 2002] in h can influence 
the results of the cloud geometrical thickness determina- 
tion in a great extent. To illustrate this, we plot the ratio of 
derivatives R, to R; (or Y = R,/R) in Figure 8. The 
parameter Y gives the relationship of the error in the 
determination of the cloud geometrical thickness A/ and 
the error in the cloud top height Ah. Namely, we have: 
Al = —YAh. 

[64] We see that the error of the cloud geometrical 
thickness determination is 1.5—3.0 times larger than the 
error with which the cloud top height is known (Y € [—3.0, 
—1.5]; see Figure 8). It means that if the cloud top height is 
obtained with the accuracy, e.g., 1.0 km, then the cloud 
geometrical thickness cannot be obtained with the accuracy 
better than 1.5—3.5 km, which is clearly too large for most 
of applications. This emphasizes the problems that can be 
accounted, e.g., in the interpretation of the GLI data with 
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on the cloud top height at the solar zenith angle 


60 degrees and the nadir observation for various values of / and 7T. Results for resolutions correspondent 
to MERIS and GLI instruments are shown (see text). 
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Figure 9. The same as in Figure 8 but for various solar angles at / = 0.5 km and t = 30. 


respect to the cloud top height determination [Nakajima et 
al., 1998]. 

[6s] Note that the cloud geometrical thickness can be 
easily obtained from the reflectance measurements in the 
oxygen A band if the value of cloud top height is known 
with a high precision (e.g., from space lidar measurements 
[Poole et al., 2003]). 

[66] The dependence of the derivative dl/dh on the solar 
angle is given in Figure 9. It shows that the high Sun cases 
are preferable for the remote sensing of the cloud geomet- 
rical thickness (for only approximately known value of the 
cloud top height). 

[67] The problem of the determination of both parameters 
can be partially solved if one uses the neural network 
approach and trains the inverse algorithm for typical cloudy 
conditions [Fischer et al., 2000]. However, this also can 
lead to large retrieval errors if these conditions are not met 
for a given experiment. The main problem here is a low 
spectral resolution both GLI and MERIS (see Figure 2), 
which makes the solution of the problem highly difficult. 

[68] This once more emphasizes that a simple Lambertian 
reflector model cannot be applied for accurate cloud top 
height retrievals. The value of / is an important parameter, 
which cannot be neglected. 

[69] The difficulty, however, can be avoided, if one uses 
the measurements from GOME or SCIAMACHY, which 
have much higher spectral resolution that those of GLI and 
MERIS (see Table 1). These advanced instruments are also 
characterized by a larger number of spectral points mea- 
sured in different regions of the oxygen A absorption band. 
The shortcomings of GOME or SCIAMACHY as far as 
cloud research is concerned are related to their large foot- 
prints. Therefore instruments with high spatial resolution 
and also with high spectral resolution in the oxygen A band 
can improve substantially our knowledge of the terrestrial 
cloudiness. Note that such an instrument (the Orbiting 
Carbon Observatory (OCO)) has already been designed 
[Kuang et al., 2002]. The OCO should be launched in 2007. 

[70] Clearly the results presented above can be obtained 
also for GOME and SCIAMACHY by a special choice of 
channels for the retrieval (correspondent to channels of GLI 


or MERIS, see Figure 2). There is even possibility of the 
optimization of channels for the retrieval to account for 
the spectral behavior of the error of our model. Note that the 
error of the retrieval of A for a known / for GOME or 
SCIAMACHY instrument using measurements for all 
wavelengths in the oxygen A band is close to data given 
for GLI. This is because in both cases almost all band is 
used. In principle, the error can be optimized with respect to 
the wavelength selection. We do not consider this possibility 
here, however. 

[71] The superiority of GOME and SCIAMACHY instru- 
ments is due to the possibility to infer accurately both 
parameters (A, /) from measurements in the oxygen A band. 
This is because instead of 1 measurement point in the O, A 
band, we have approximately 66 (for GOME) and 44 (for 
SCIAMACHY) spectral points (see above). 

[72] There are several possibilities for the development of 
the retrieval algorithm in this case. We chose here the 
simplest one, which is based on the technique, described 
in the previous section. Namely, we assume a small geo- 
metrical thickness of the cloud equal, e.g., to 100 m and 
perform the retrieval to find the value of the cloud top 
height A as described in section 3.1 using all spectral points 
available. Then the value of the maximal difference 6, 
between the calculated O) A band spectrum and measured 
spectrum is derived, taking into account all spectral points. 
The value of / is increased and the procedure is repeated till 
the minimum of the function 6; (/) is found. 

[73] Therefore the solution of the two-parameter problem 
is reduced to the sequence of solutions of the one-parameter 
problem of the finding / at given value of J, (k= 1, ... N) 
(see section 3.1). Namely it follows that 


The iteration process used to find h (/,) allows to reduce 
the influence of the linearization errors [Rodgers, 2000] in 
the case when the initial approximation h© (/,) is far away 
from the solution of the problem. Note that the process of 
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Figure 10. The illustration of the performance of simultaneous cloud top height and cloud geometrical 
thickness semianalytical retrieval algorithm. Measurements were simulated using SCIATRAN [Rozanov 
et al., 2002]. The random error of 0.3% was added to measurements. Shaded areas show cloud 
parameters for a simulated case (“true clouds”). Dashed lines show retrieved clouds. The number 1 is 
related to the retrieval of cloud top height for a known value of /. The number 2 is related to the case, 
when both A and / are not known. Calculations are given for the GOME resolution at T = 30 and nadir 
observation. The solar zenith angle is equal to 53 degrees. Different panels correspond to different “true 


clouds”’ positions. 


iterations is terminated when the difference between 
retrieved values of cloud top heights becomes smaller than 
50 m. Then we calculate the value of the error: 


stay = pa) ~ (0.4) 


for a spectral range studied. 

[74] We emphasize that the algorithm described was 
developed to demonstrate the dependence of the retrieved 
cloud top height on the cloud geometrical thickness. Other 
possible solutions of the inverse problem can be used 
[Rodgers, 2000]. 

[75] Let us consider now Figure 10, where we assumed a 
random error of 0.3% for simulated radiances. Note that 
these radiances are also in error (typically below 5%) owing 
to our model inaccuracy as compared to SCIATRAN data. 
The vector R is given by a column of numbers for the 
GOME as specified above. 

[76] Note that the case 1 in Figure 10 gives the retrieved 
cloud top height if the value of / is exactly known. The case 


(24) 


2 corresponds to the case when both / and / are not known 
and the pair (h, /) is chosen using the location of the 
minimum of the function 6; (Z). 

[77] Data given in Figure 10 allow to make the following 
conclusions: 

[78] 1. The iteration algorithm (see equation (23)) finds 
the cloud top height / with the accuracy better than 100 m if 
the cloud geometrical thickness / is known. 

[79] 2. The error of the cloud top height determination 
increases up to 250 m if both h and l are unknown 
parameters of the problem. The error of the cloud geomet- 
rical thickness is in the range 200—600 m, depending on a 
particular case studied. Note that the influence of random 
errors is very profound in this case (even for their small 
values). This is because this error displaces the minimum of 
the function 6, (/,) considerably (see equation (24)). 

[s0] 3. The dependence of the retrieved value of A on / can 
be approximated by a linear function. Our estimations lead 
to values of di/dh in the range 1.6—1.7. This is in an 
accordance with results shown in Figures 8 and 9. 
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LWC profiles (see text for explanation). 


[s1] Summing up, we see that the retrieval accuracy is 
quite high if synthetic data are used. In particular, the value 
of h is found with an accuracy better than 250 m if both 
h and / are unknown parameters of the problem. The 
important question is the performance of the algorithm for 
the case of satellite measurements. Those are considered in 
a separate publication [Rozanov et al., 2004]. 


4. Errors of Retrievals as Applied to Vertically 
Inhomogeneous and Multilayered Clouds 


[s2] The previous section was devoted to the determi- 
nation of the parameters h, l, assuming that a cloud field 
under study can be represented by a single homogeneous 
layer. Such an assumption is never valid. The concentra- 
tion and size of droplets in clouds change with height 
[Feigelson, 1984]. Also, multilayered cloud systems are 
more common than a simple case of a single cloud layer 
[Grechko et al., 1973, 1975; Platnick et al., 2003]. Many 
clouds contain also ice crystals of various shapes and not 
only spherical water droplets. This certainly leads to 
biases of retrieved parameters if a single liquid water 
cloud with the constant liquid water content (LWC) is 
assumed in the retrieval procedure. A brief study of these 
biases using several numerical experiments is given in this 
section. 


4.1. Influence of the Vertical Profile of the 
Liquid Water Content 


[s23] We choose four liquid water content profiles for the 
study of the influence of these profiles on the retrieval 
procedure. They are given in Figure 11. The dimensionless 
height is defined as v = (h — z)(h — b), where z is the 
vertical coordinate, changing from z = b at the cloud bottom 
to z = hat the cloud top. The profile scaling factor s (v) (see 


1.0 
Profile—scaling factor 


1.5 2.0 


The dependence of the profile-scaling factor on the dimensionless height for four types of 


Figure 11) determines the cloud optical thickness T and 
liquid water path W with following equations: 


h 


T= (ext) | s(z)dz, (25) 
b 
h 
W = (LWC) | s(z)dz, (26) 
| 


where (Oe) and (LWC) are the average values of the 
extinction coefficient and LWC inside the cloud. Therefore 
the liquid water profile inside the cloud is determined as the 
product s (z) (LWC). 

[s4] The profile of the type 1 (see Figure 11) is found in 
most of cloud layers [Feigelson, 1984]. It gives maximum 
of the LWC near the top of the cloud. The profile of the 
second type has been used in calculations of the influence of 
clouds on climate [Marchuk et al., 1986]. The profile of the 
third type is the superposition of profiles 1 and 2. The 
profile 4 corresponds to the vertically homogeneous layer 
used above. Note that we neglect the possible variability of 
cloud droplets sizes with height here. Also, it is assumed 
that ice crystals are absent. 

[ss] Let us consider now the results of the retrieval of 
parameters h and / with our algorithm, assuming different 
types of liquid water content profiles for synthetic data. The 
inverse algorithm, however, is the same as for the case of 
the LWC = const. It means that we consider the case when 
the forward model does not match the realistic cloud 
system. 

[s6] For this, we have created synthetic top of atmo- 
sphere reflectance data Rmes using SCIATRAN [Rozanov 
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(a) The influence of the vertical profile of the liquid water content (see Figure 11) on the 


accuracy of the determination of values h, /. The spectral resolution corresponds to GOME. It is assumed 
that cloud optical thickness is equal to 30 and the solar zenith angle is 53 degrees. The random error is 
equal to 0.3%. The observation is in the nadir direction. Big squares give the position of the “true” cloud 
in the coordinates (A, /). Different symbols give results of the retrieval under different choice of the LWC 
in a cloud (see Figure 11). Numbers in Figure 12a correspond to numbers in Figure 11. The cases 1, 2, 
and 3 correspond to clouds with values of (A, /) equal to (3.0, 0.5), (4.0, 1.0), and (3.0, 2.0) km, 
respectively. (b) Same as in Figure 12a, but the synthetic data have been calculated using the 
semianalytical technique described above (and not SCIATRAN as in Figure 12a). 


et al., 2002] for the vertical profile of the type 1. The 
retrieval results using the synthetic spectra are summarized 
in Figure 12 for three cases, having the parameters (h, /) 
equal to (3.0, 0.5) (case 1), (4.0, 1.0) (case 2), and (3.0, 
2.0) (case 3). The random error € = 0.3% is added to Rines- 
It follows from Figure 12a that the error of the cloud top 
height determination increases up to 600 m if it is assumed 
that LWC = const. The difference of Figure 12b from 
Figure 12a is that the value of Rmes was calculated not 
using the SCIATRAN but the asymptotic theory given 
above. This allows to model the influence of the random 
error and a wrong LWC profile on results of retrievals. 
The errors of the radiative transfer model are excluded in 
this case. 

[s7] It follows that the influence of the LWC profile is 
more important for clouds having larger geometrical thick- 
ness. For instance, the use of the condition LWC = const 
in the inversion procedure gives the error of 400 m for the 
case 1 (J = 0.5 km). The error is equal to 750 m for the 
case 3 (J = 2 km). The comparison of Figures 12a and 12b 
shows that the use of the semianalytical cloud retrieval 
scheme does not lead to the considerable increase of the 
retrieval error. 

[ss] Summing up, the use of the assumption LWC = const 
in the retrieval procedure for the case of realistic vertically 
inhomogeneous clouds can lead to potentially large errors 
(up to 600—800m in h and up to 1 km in J). Clearly the 
assumption of a typical LWC profile (and not LWC = const) 


in the retrieval procedure for a single cloud case will lead to 
the decrease of errors mentioned above. 


4.2. Influence of the Second Cloud Layer 


[s9] Let us consider now the case of two-layered cloud 
systems and their influence on the retrieval procedure on the 
basis of the assumption of a homogeneous single cloud 
layer. Generally, the presence of multilayered cloud systems 
[Grechko et al., 1973, 1975] leads to a great variety of 
different situations, which differ not only by the number, 
position, geometrical and optical thickness of cloud layers, 
but also by a possibility to have water droplets or exclu- 
sively ice crystals in different cloud layers. 

[90] Let us consider now the results of the retrieval of 
parameters h and / with our algorithm, assuming different 
types of two-layered cloud systems. Namely, we will 
consider the case of two liquid water cloud layers and 
separately the case of an ice cloud above liquid one. Such 
situation is a typical one for terrestrial atmosphere. 

[91] We start from the case of liquid clouds. For this, we 
have created synthetic top of atmosphere reflectance data 
Rimes using SCIATRAN [Rozanov et al., 2002] for the two- 
layered liquid cloud system. Also, we have assumed that the 
vertical profile of the LWC in both clouds correspond to the 
type 1, given in Figure 11. The results of retrievals assum- 
ing the single-layer cloud system with LWC = const are 
given in Figures 13 and 14 for various values of optical 
thickness of upper and lower cloud. The analysis of data 
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given in these figures, which differ only by the distance Z 
between two cloud systems (Z = 600 m for Figure 13 and 
Z = 2 km for Figure 14), allows us to reach following 
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Figure 13. The results of the retrieval of A and / of the cloud system composed of two water cloud 
layers. The distance between layers Z is equal to 600 m. Big squares give the position of the “true” cloud 
in the coordinates (h, /). Different symbols give results of the retrieval under different choices of the 
LWC in a cloud (see Figure 11). Numbers under curves give the cloud optical thickness of the lower and 


upper cloud layer, respectively. 


The solid line (A = 1) gives the position of the cloud having the value of h 


equal to that of /. Dark area corresponds to the upper boundary of the upper cloud +1 km. 


conclusions: 
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Figure 14. The same as in Figure 13 but Z = 2 km. 
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[92] 1. The assumption LWC = const and a single cloud 
layer in the retrieval procedure gives the cloud top height of 
the vertically inhomogeneous two-layered cloud system 
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Figure 15. The results of the retrieval of A, Z of the cloud system, composed of the upper ice cloud and 
lower water cloud. The distance between clouds Z is equal to 600 m. Other symbols have the same 


meaning as in Figure 13. 


[93] 2. The geometrical thickness of the retrieved 
“effective” single cloud depends on the geometrical thick- 
ness of the layer between two clouds. The geometrical 
thickness of the retrieved “effective” cloud is equal to 
4.3—4.8 km at Z = 600 m for considered cases (see 
Figure 13). This is close to the geometrical thickness of 
the layer between the cloud system top height and under- 
lying surface. Interestingly, the value of / increases up to 
6-7 km at Z = 2 km (see Figure 14). Then we have the 
unphysical condition / > h. This can be explained as 
follows. Namely, the cloud geometrical thickness deter- 
mines the amount of light absorption in the cloud in the 
oxygen A band. Therefore the increase of Z leads to the 
increase of the absorption in the two-layered cloud system 
owing to the increase of amount of oxygen between cloud 
layers. This leads to the values of / larger than A. This 
unphysical result is due to the use of a single cloud with 
LWC = const in the retrieval procedure. In fact, the two- 
layered system should be used in the retrieval procedure. 
However, it is difficult to decide if a single cloud layer or 
a multilayered cloud system is responsible for a given 
TOA reflectance. By the way, if the retrieval gives / > h, 
then it is most probably that we have a multilayered case. 
Therefore our method can provide an indirect information 
on cloud system vertical structure. This by itself is a 
valuable information. Further research is needed in this 
direction, however. 

[04] 3. The use of other types of the LWC profiles in the 
retrieval procedure for the two-layered case leads to sys- 
tematically lower estimates of h and / as compared to 
results, obtained under assumption LWC = const. The 
difference reaches 1 and 2 km at Z = 600 m for h and J, 
respectively. It is even larger (%2.5 km) at Z = 2 km (see 
Figure 14). 


[95] 4. The optical thickness T* of the two-layered 
system, retrieved using approach giving by Kokhanovsky 
et al. [2003], gives the value approximately equal to the 
sum of optical thicknesses of both layers. 

[96] In conclusion, let us consider the results of the 
retrieval of parameters / and / with our algorithm, assuming 
the case of an ice cloud above liquid one. For this, we have 
created synthetic top of atmosphere reflectance data Rimes 
using SCIATRAN [Rozanov et al., 2002] for this two-layered 
cloud system. The phase function of the ice cloud was 
modeled using the hexagonal prisms in according with the 
model of Mishchenko et al. [1999]. The height dependence of 
the phase function was neglected. Also, we have assumed 
that the vertical profile of the liquid/ice water content in both 
clouds correspond to the type 1, given in Figure 11. 

[97] The results of retrievals assuming the single-layered 
liquid cloud system with LWC = const are given in 
Figures 15 and 16 at Z = 600 m and Z = 3.6 km, 
respectively. The optical thickness of the upper cloud layer 
was assumed to be between 1 and 10. Figures 15 and 16 are 
similar in many ways to Figures 13 and 14. This indicates 
that the thermodynamic state of the upper cloud influence 
the retrieval insignificantly. There are some differences, 
however. First of all, the crystalline clouds are usually 
optically thin. This leads to the increase of the retrieval 
error. In particular, the error is equal to ~1.3 km at Z = 
600 m and ~2.5 km at Z = 3.6 km assuming a single cloud 
layer with LWC = const in the retrieval procedure. 

[98] Second, the optical thickness of the retrieved cloud 
system T is approximately two times larger than the sum of 
optical thicknesses of two separate cloud layers T* in this 
case. Remember that we had T ~ t* in the case of the liquid 
cloud system. The reason for this is that the existence of ice 
phase in the cloud is ignored in the retrieval procedure. 
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Figure 16. The same as in Figure 15, except at Z = 3.6 km. 


Therefore the information on the cloud thermodynamic state 
is an essential for a correct run of the retrieval code. 


5. Conclusion 


[99] We presented here a simple and flexible asymptotic 
algorithm for the determination of cloud top height from 
measurements in the oxygen A band. Our technique can be 
also applied to other gaseous absorption bands. It relied on 
the semianalytical solution of the radiative transfer equation, 
which is valid for clouds having the cloud optical thickness 
larger than 5 [Kokhanovsky and Rozanov, 2003, 2004]. Note 
that thick clouds with T > 5 occur quite often [Trishchenko 
et al., 2001]. The comparison of our algorithm with 
other retrieval methods is given in a separate publication 
[Rozanov et al., 2004]. 

[100] The accuracy of the retrieval algorithm for the cloud 
top height determination if a single water cloud with LWC = 
const and known geometrical thickness is present in the 
scene under study is better than ~200 m. This error 
increases by ~50 m if the value of the cloud geometrical 
thickness is also unknown parameter. Then / can be also 
estimated. However, the error of the cloud geometrical 
thickness determination is generally larger (up to 3 times; 
see Figures 8 and 9) than that for the case of the cloud 
altitude determination. 

[101] We also study the performance of our algorithm if a 
two-layered cloud system or a cloud with LWC ¥ const is 
present in the scene under study. Then errors of the retrieval 
increase considerably if a single cloud layer with the LWC = 
const is assumed in the retrieval procedure. Therefore it is of 
a great importance to select a right cloud model in the 
retrieval algorithm. This is a difficult problem, however, for 
measurements in real atmospheric conditions. 

[102] Note that values of the cloud optical thickness, 
effective radius of droplets, and liquid water path can be 
found outside of the oxygen bands, using recently developed 


cloud retrieval algorithms [Kokhanovsky et al., 2003; 
Platnick et al., 2003], which are based on the similar 
assumptions as in the asymptotic algorithm described above. 

[103] All studies presented above have been performed 
for the black underlying surface. The effects of polarization 
and non-plane-parallel geometry are neglected. The uncer- 
tainty of the ground reflectance (especially over land) can 
increase the retrieval error, so surface albedo databases 
should be used to avoid the problem. Our estimations (not 
presented here) show that the uncertainty in albedo of 50% 
(e.g., change of surface albedo from 0.2 to 0.3) results in 
150 m (or smaller) error in the cloud top height determina- 
tion if the cloud optical thickness is larger than 20. For 
thinner clouds the influence of incorrect albedo increases 
the retrieval error substantially. This confirms that accurate 
estimations of the surface albedo over land should be used 
as an input to our cloud retrieval algorithm. This problem is 
only of a minor importance for measurements over ocean. 

[104] The combination of the cloud retrieval algorithm 
mentioned above [Kokhanovsky et al., 2003] with that given 
in this paper allows to obtain a number of other cloud 
characteristics of interest (e.g., the liquid water content 
and cloud extinction coefficient) and constitutes a newly 
developed Semi-Analytical Cloud Retrieval Algorithm 
(SACURA). Note that the SACURA is much more faster 
and flexible than approaches based on the exact solution of 
the radiative transfer equation. 


Appendix A: Auxiliary Functions 


[105] Here we introduce without derivations auxiliary 
functions, which have been used in the main text of the 
paper (see Table 2). Details are given by Kokhanovsky and 
Rozanov [2003]. 

[106] Various parameters in Table 2 have the meaning 
described below [see also Kokhanovsky and Rozanov, 
2003]. The value of R% (€, n, ~) gives us the reflection 
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Table 2. Auxiliary Functions [Kokhanovsky and Rozanov 2004; Kokhanovsky et al., 2003] 


Auxiliary Function Formula 
Re REO, Vo, pJexp (—(1-0.05y)u(9, Do, P)) —i* exp (~x —y)Ko (Oo)Ko (9) 
u Ko(0)Ko(Yo) 
R? (0,009) 
Ko (9) 31 +2 cos Ù) 
q te — (4.86— 13.08 cos Jy cos 0 + 12.76cos” By cos” O)exp O/T 
sinhy 
te saher OD], 
fo exp (—y (1 — 0.05y)) — ¢* exp (—x —y) 
1 1 
Ya wos EdE f ndnF(& np (& n) 
0 0 
1—exp[-r4(€1 +1°1)] 
F F etn 
Ra gia f ElOsca OP” ©) + o Op E, exp [TEE + 
h 


function of a semi-infinite nonabsorbing cloud. It has the 
following simple form for close to nadir observations 
[Kokhanovsky, 2004]: 


A+B(E +n) + Cn +f) 
4(€ +n) 


RY (& 1,9) = 


where f (0) = p (0) — P (& 1), 9 = arccos (—€n + 
y (L= ©) (I — W) cos 9), € = cos Yo, n = loos 0], 2 &, n) 
= fp Ode, and A ~ 3.944, B ~ —2.5, C ~ 10.664. Here 


0, 
p(®) is the cloud phase function. The accuracy of this 
approximation is studied in detail by Kokhanovsky [2004]. 


1—wo 


[107] Parameters x = yd, y= 4 , [3g 4 = 3 (l-8) Te 
are needed to describe the reflection function of a cloud at 
absorbing wavelengths. Here g is the asymmetry parameter 
[Kokhanovsky, 2001], wo is the single scattering albedo, and 
To is the cloud optical thickness. 

[108] We found approximating the function r, (T^) at the 
asymmetry parameter g = 0.77 and 7^ < 0.3: 


ra(T) = 0.26574 — 0.865 (14)? 41.164 (74)°, 


where 7^ is the aerosol optical thickness beneath the cloud. 

[109] Other parameters in Table 2 are defined as follows: 
(1) o4,(z) and o£, (Z) are total aerosol and Rayleigh 
scattering coefficients, respectively; (2) p“ (0, z) and p° 
(8) are aerosol and Rayleigh phase functions, respectively; 
(3) T (Z) is the optical depth of the atmosphere above z along 
the vertical axis OZ (see Figure 1), which includes the 
contribution from both molecular and aerosol scattering and 
absorption; and (4) H is the top of atmosphere height, 
assumed to be equal to 60 km in this study. 

[110] The single scattering albedo wọ is defined by the 
following equation: 


Oabs 
Wo = 1- 


? 
Oext 


where Oabs — Tás at Gite T Ths and Oext 7 Tka + ee F ELi 
where indices A, G, and C show aerosol, gas, or cloud 
contribution, respectively, to extinction Oext or absorption 
Oabs Coefficients. 


[111] Let us assume that outs = of, = 0. Then we have 


G 

fo) 
wy = 1 — abs, 

Oext 


Thus the value of wọ changes with the height both inside 
and outside of a cloud. Formula for the cloud reflection in 
Table 2, however, is applicable only for the case of a 
vertically homogeneous layer. The dependence wọ (z) is not 
particularly strong in the area, where cloud is present. 
Therefore we adopt here the model of an “effective 
homogeneous layer” [Chamberlian, 1965; Yanovitskij, 
1997]. In this case one should find the height inside the 
cloud at which the value of wo should be taken for 
calculations. Details of this are given by Kokhanovsky and 
Rozanov [2003]. Note that we also have accounted for the 
variation of the cloud liquid water content with the height. 
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